/**
  ETFAna project, Anyang Normal University && IMP-CAS
  \class ETFAuxAna
  \brief a worker class auxiliary to ETFAnalyzer, for refined ana in a global scope
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2024-06-23
  \date 2024-06-23 last modified
  \attention
  changelog
  <table>
  <tr>  <th>Date         <th>Author      <th>Description                   </tr>
  <tr>  <td>2024-06-23   <td>Asia Sun    <td>file created                  </tr>
  </table>

  \copyright Copyright (c) 2021-2024 Anyang Normal U. && IMP-CAS with LGPLv3 LICENSE
*/

#include <TTree.h>
#include "ETFAuxAna.h"
#include "ETFDetectorManager.h"
#include "ETFTOFWall.h"
#include "ETFDCTa.h"
#include "ETFDC.h"
#include "ETFMUSIC.h"
#include "ETFT0Array.h"
#include "ETFTimer2.h"
#include "ETFTimer4.h"
#include "ETFHPTDCCh.h"
#include "ETFMUSIC.h"
#include "ETFMath.h"
#include "ETFAmpCh.h"
#include "ETFPlaHit.h"
#include "ETFDCHit.h"
#include "ETFVPID.h"
#include "ETFRun.h"
#include "ETFMsg.h"

ETFAuxAna::ETFAuxAna(){}

ETFAuxAna::~ETFAuxAna(){}

void ETFAuxAna::Init(){
  auto dm = ETFDetectorManager::Instance();
  fT1p = dynamic_cast<ETFTimer2 *>(dm->GetDetector("tstop"));
  fdcTa[0][0] = dynamic_cast<ETFDCTa *>(dm->GetDetector("dcTaU0"));
  fdcTa[0][1] = dynamic_cast<ETFDCTa *>(dm->GetDetector("dcTaU1"));
  fdcTa[1][0] = dynamic_cast<ETFDCTa *>(dm->GetDetector("dcTaD0"));
  fdcTa[1][1] = dynamic_cast<ETFDCTa *>(dm->GetDetector("dcTaD1"));
  fdc[0] = dynamic_cast<ETFDC *>(dm->GetDetector("dc0"));
  fdc[1] = dynamic_cast<ETFDC *>(dm->GetDetector("dc1"));
  fdc[2] = dynamic_cast<ETFDC *>(dm->GetDetector("dc2"));
  fT1v = dynamic_cast<ETFTimer2 *>(dm->GetDetector("tstop0_v"));
  fVETO = dynamic_cast<ETFTimer4 *>(dm->GetDetector("veto_v"));
  fT0p = dynamic_cast<ETFTimer2 *>(dm->GetDetector("tstart2"));
  fT0v = dynamic_cast<ETFTimer2 *>(dm->GetDetector("tstart2_v"));
  fTOFW = dynamic_cast<ETFTOFWall *>(dm->GetDetector("tofwB"));
  if(!fTOFW) fTOFW = dynamic_cast<ETFTOFWall *>(dm->GetDetector("tofw"));
  fT0F1 = dynamic_cast<ETFTimer4 *>(dm->GetDetector("tstart4"));
  fT0Arr = dynamic_cast<ETFT0Array *>(dm->GetDetector("t0arr"));
  fMU[0] = dynamic_cast<ETFMUSIC *>(dm->GetDetector("mu0"));
  fMU[1] = dynamic_cast<ETFMUSIC *>(dm->GetDetector("mu1"));
} // end member function Init

void ETFAuxAna::Reset(){
  atof1 = -999.;
  aaoz0 = -999.;
  az0 = abeta0 = abg0 = -1.;
} // end member function Reset

void ETFAuxAna::Register(){
  TTree *tree = ETFRun::Instance()->GetTree();
  tree->Branch("atof1", &atof1, "atof1/D");
  tree->Branch("abeta0", &abeta0, "abeta0/D");
  tree->Branch("abg0", &abg0, "abg0/D");
  tree->Branch("az0", &az0, "az0/D");
  tree->Branch("aaoz", &aaoz0, "aaoz0/D");
} // end member function Register

void ETFAuxAna::Ana(){
  Reset();

  //////////////////////// SECONDARY BEAM PID ////////////////////////////////
  // get tof1 //
  if(fT0p && fT1p && fT0p->fired() == 4 && fT1p->fired() == 4)
    atof1 = fT1p->time() - fT0p->time();
  // get DeltaE0, and hit pos at ETF and F1 //
  double mu0 = -999., dcTaU1_wireHitPos[2] = {-999., -999.}, xf1 = -999.;
  short mu0v[10]{};
  if(fMU[0]) mu0 = fMU[0]->dE(); // MUSIC0
  for(int i = 0; i < fMU[0]->nunit(); i++) mu0v[i] = fMU[0]->dE(i);
  if(fdcTa[0][1]){
    dcTaU1_wireHitPos[0] = fdcTa[0][1]->WireHitPosX(); // hit X in DCTaU1
    dcTaU1_wireHitPos[1] = fdcTa[0][1]->WireHitPosY(); // hit Y in DCTaU1
  } // end if
  if(fT0Arr) xf1 = fT0Arr->GetOptimalHitX(); // hit X in T0 Array @F1
  // get beta*gamma //
  // static const double BRHO0 = 4.9896745, BRHOC = BRHO0 * ETFMath::CC(); // 1.6058799
  static const double LL = ETFVPID::L / ETFMath::c0();
  abeta0 = LL / atof1;
  if(abeta0 > 0. && abeta0 < 1.) abg0 = ETFMath::BetaGamma(abeta0);



  ///--- now we have the ingredient to calibrate z0 ---///
  if(mu0 > 0. && dcTaU1_wireHitPos[1] != -999. && atof1 > 0.){
    // > jun2024 - normal runs
    // az0 = sqrt(mu0*exp(0.0013*(dcTaU1_wireHitPos[1]+30.)) / (0.0027*pow(atof1-121.,2)+7.))*53.9046-0.480279;
    // > jun2024 - the low z run
    az0 = sqrt((mu0v[2]*exp(0.0015*(dcTaU1_wireHitPos[1]+30.))+mu0v[3]*exp(0.0022*(dcTaU1_wireHitPos[1]+30.)))/2000./(0.0016*pow(atof1-100., 2.)+5.))*50.7183-1.34103;
  } // end if
  ///--- now we have the ingredient to calibrate aoz0 ---///
  if(-999. != xf1 && -999 != dcTaU1_wireHitPos[0] && abg0 > 0.){
    // > jun2024 - normal runs
    // aaoz0 = BRHOC/bg0*(1.+0.000194047*hitX-0.0011855*x_f1)+0.19894667;
    // > jun2024 - the low z run
    aaoz0 = 1.6058799/ETFMath::BetaGamma(86.33/(atof1+13))*(1.+0.0003*dcTaU1_wireHitPos[0]-0.00105*xf1-0.000002*xf1*xf1)*1.25631-0.249302;
  } // end if
  //////////////////////// SECONDARY BEAM PID ////////////////////////////////

  /////////////////////////// OTHER ANA //////////////////////////////////////
  /////////////////////////// OTHER ANA //////////////////////////////////////
} // end member function Ana
